!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Carlo de Falco                   <carlo.defalco@gmail.com>

module mod_interp1d

!General comments: provides some basic routines for 1d interpolation.
!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------

 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &

   mod_interp1d_constructor, &
   mod_interp1d_destructor,  &
   mod_interp1d_initialized, &
   ! lookup values in an ordered list
   lookup,                   &
   ! interpolate linearly in 1d 
   linterp1d

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 ! private members

! Module variables

 ! public members
 interface linterp1d
    module procedure scalar_linterp1d, vect_linterp1d
 end interface linterp1d

 interface lookup
    module procedure scalar_lookup, vect_lookup
 end interface lookup

 logical, protected ::               &
   mod_interp1d_initialized = .false.

 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_interp1d'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

  subroutine mod_interp1d_constructor ()
    character(len=*), parameter :: &
         this_sub_name = 'constructor'
    
    !Consistency checks ---------------------------
    if ((mod_messages_initialized.eqv..false.) &
         .or. (mod_kinds_initialized.eqv..false.)) then
       call error (this_sub_name,this_mod_name, &
            'Not all the required modules are initialized.')
    endif
    if (mod_interp1d_initialized.eqv..true.) then
       call warning (this_sub_name,this_mod_name, &
            'Module is already initialized.')
    endif
    !----------------------------------------------
    
    mod_interp1d_initialized = .true.
  end subroutine mod_interp1d_constructor
  
  !-----------------------------------------------------------------------
  
  subroutine mod_interp1d_destructor ()
    character(len=*), parameter :: &
         this_sub_name = 'destructor'
    
    !Consistency checks ---------------------------
    if (mod_interp1d_initialized.eqv..false.) then
       call error (this_sub_name,this_mod_name, &
            'This module is not initialized.')
    endif
    !----------------------------------------------

    mod_interp1d_initialized = .false.
  end subroutine mod_interp1d_destructor
  
  !-----------------------------------------------------------------------
  pure subroutine vect_lookup (array, values, indices)
    ! Given an ordered array "array" and a list of values "values"
    ! find the index of the smallest element in array that is strictly 
    ! less than each element of "values"
    
    real(wp), intent(in)  :: array(:)    ! ordered list of real values
    real(wp), intent(in)  :: values(:)   ! values to lookup
    integer,  intent(out) :: indices(:)  ! indices of values in list
    integer :: m, n, i
    
    m = size(array)
    n = size(values)
    
    do i = 1, n
       call scalar_lookup (array, values(i), indices(i))
    end do
    

  end subroutine vect_lookup
  
  !-----------------------------------------------------------------------
  pure subroutine scalar_lookup (array, val, idx)
    ! Given an ordered array "array" and a value "val"
    ! find the index "idx" of the smallest element in array that is strictly 
    ! less than  "val"
    
    real(wp), intent(in)  :: array(:)    ! ordered list of real values
    real(wp), intent(in)  :: val         ! valu to lookup
    integer,  intent(out) :: idx         ! index of value in list
    integer :: m
    
    m = size(array)
    idx = 0;
    main_loop: do while (idx .lt. m)
       if (val .ge. array(idx + 1)) then
          idx = idx + 1
       else
          exit main_loop
       end if
    end do main_loop
  end subroutine scalar_lookup

  !-----------------------------------------------------------------------
  
  pure subroutine vect_linterp1d (xin, yin, xout, yout, indices)
    ! Interpolate the piecewise linear function described by the xin, yin
    ! at the nodes xout. if a lookup
    ! of the subintervals where the evaluation points xout are located has 
    ! already been performed, then pass the resulting indices to save time 

    real(wp), intent(in)  :: xin(:), yin(:), xout(:) ! ordered list of real values
    real(wp), intent(out) :: yout(:)                 ! indices of values in list
    integer,  intent(in), optional :: indices(:)     ! output of lookup (xin, xout, indices)
    integer :: idx(size(xout))
    integer :: ii, m, n

    m = size (xin)
    n = size (xout)

    if (.not. present (indices)) then
       call lookup (xin, xout, idx)
    else
       idx = indices
    end if

    do ii = 1, n
       call scalar_linterp1d (xin, yin, xout(ii), yout(ii), idx(ii))
    end do
    
  end subroutine vect_linterp1d

  !-----------------------------------------------------------------------

  pure subroutine scalar_linterp1d (xin, yin, xout, yout, idx)
    ! Interpolate the piecewise linear function described by the xin, yin
    ! at the single node xout. if a lookup
    ! of the subintervals where the evaluation points xout are located has 
    ! already been performed, then pass the resulting index idx to save time 

    real(wp), intent(in)  :: xin(:), yin(:), xout    ! ordered lists of real values
    real(wp), intent(out) :: yout                    ! interpolated values
    integer,  intent(in), optional :: idx            ! output of lookup (xin, xout, indices)
    integer :: idxi
    integer :: m

    m = size (xin)

    if (.not. present (idx)) then
       call scalar_lookup (xin, xout, idxi)
    else
       idxi = idx
    end if

    if (idxi .lt. 1) then
       yout = yin(1)
    elseif (idxi .ge. m) then
       yout = yin(m)
    else
       yout = &
            (yin(idxi) * (xout - xin(idxi+1)) - &
            yin(idxi+1) * (xout - xin(idxi))) / &
            (xin(idxi) - xin(idxi+1))
    end if

  end subroutine scalar_linterp1d
  !-----------------------------------------------------------------------
  

end module mod_interp1d

